退化过程模型综述

庄亮亮、吴温慧

近年来，基于退化的可靠性技术在模型、方法和应用等方面得到快速发展。其中，在基于退化的可靠性模型方面，以随机过程理论为基础，根据工程需要，并考虑模型的简明性、实用性和适用性，已经提出了多种类型的模型，包括：

基于 Wiener 过程的模型

基于 Gamma 过程的模型

基于逆高斯过程的模型

在基于退化的可靠性建模方法，以数理统计理论为基础，针对模型和退化数据的类型，研究矩估计、极大似估计及 EM 算法、Bayes 估计及 MCMC算法、基于滤波的状态估计方法，以及非参数和半参数方法等，解决模型辨识和修正问题;

采用似然比检验、Bayes 因子分析等方法进行模型验证;

采用均方误差（MSE）、Akaike 信息准则(AIC)、Bayes 信息准则(BIC)、偏差信息准则(DIC)等进行模型优良性检验。

在基于退化的可靠性技术应用方面，针对机械零部件（如轴承、润滑系统）、半导体器件（如功率MOS器件）、机电部件（如感应电动机）、光电器件(如太阳电池、激光器、LED)、电子器件（如电容、蓄电池）等的退化失效过程，开展退化过程建模、可靠性评估、剩余寿命预测、可靠性试验设计：特别是加速退化试验方案设计的研究。

Wiener 过程由于其简单的结构和比较丰富的研究成果，成为目前退化过程建模中应用最为广泛的一种模型。Wiener 退化过程的首达时间分布具有解析形式，便于对产品寿命和可靠度的分析和计算。

Wiener 退化过程可以描述产品退化过程的时间不确定性，而且比较容易处理测量数据存在误差的情况。

通过对经典 Wiener 过程的参数引入随机性，使得 Wiener 退化过程能够描述个体差异，并且一般不会给模型参数的估计带来更多的困难。

对 Wiener 退化过程及其观测过程所构成的状态空间模型，可以采用 Kalman 滤波技术进行处理，为产品在线寿命预测、视情维修决策等提供了可行的算法和实现途径。

称 \(\{X(t), t \geqslant 0\}\) 是漂移参数为 \(\mu\) 、扩散参数为 \(\sigma\) 的 (一元) Wiener 过程, 若

\(X(0)=0\)。

\(\{X(t), t \geqslant 0\}\) 有平稳独立增量。

\(X(t)\) 服从均值为 \(\mu t\)，方差为 \(\sigma^2 t\) 的正态分布。

根据上面的定义, 可以将带漂移的 Wiener 过程表示为下面的形式: \[ X(t)=\mu t+\sigma B(t) \] 式中: \(\{B(t)\}, t \geqslant 0\) 是标准 Wiener 过程或标准布朗运动过程。

根据定义, 对带漂移的 Wiener 过程, 显然有如下的性质成立:

（1）时刻 \(t \sim t+\Delta t\) 之间的增量服从正态分布, 即 \(\Delta X=X(t+\Delta t)-X(t) \sim\) \(N\left(\mu \Delta t, \sigma^2 \Delta t\right)\) 。

（2）对任意两个不相交的时间区间 \(\left[t_1, t_2\right],\left[t_3, t_4\right], t_1

基于前面所提的经典线性维纳过程，利用 rstan 进行贝叶斯分析。

\[ \begin{aligned} &\Delta X_{ij} = N(\mu \Delta t_{ij}, \sigma^2 \Delta t_{ij}),\\ &w = \frac{1}{\sigma^2} \sim Gam(a,b),\\ & \mu|w \sim N(d,\frac{c}{w}). \end{aligned} \]

其中，\(a = b = 1, d = 0, c =100\)。

准备\(\Delta t_{ij}\) 和 \(\Delta x_{ij}\)的数据，并且标志出矩阵的维数。

从性能退化过程的角度看，采用上述带漂移的 Wiener 过程模型，相当于认为同一批产品的退化过程相同，即参数心和口相同。实际也可能由于产品的原材料、生产过程等的差异，退化过程参数\(\mu\)和\(\sigma\)可能因产品个体差异而不同。为了刻画这种情景，我们给出带随机效应的 Wiener 过程模型。给出几种常用的模型：

\[ X(t)=\mu t+\sigma B(t) \]

假定 \(\mu\)为随机变量，并服从\(N(\mu_\beta,\sigma^2_\beta)\)，来描述个体退化速率的差异，\(\sigma\)是对所有个体相同的扩散参数。

文献 [5] 提出一种不同个体具有不同退化速率和扩散行为的随机效应模型，其中，认为不同单元的漂移和扩散参数都是随机变量，且服从正态-逆Gamma分布。假设

\[ \sigma^{2} \sim InvGama(\alpha, \beta), \quad \mu \mid \sigma^{2} \sim N(v, \eta \sigma^2). \]

\[ f_T(t)=\sqrt{\frac{l^2}{2 \pi t^3\left(\sigma_\beta^2 t+\sigma_B^2\right)}} \exp \left\{-\frac{\left(l-\mu_\beta t\right)^2}{2 t\left(\sigma_\beta^2 t+\sigma_B^2\right)}\right\}, t \geqslant 0 \] 对应的不可靠度函数为 \[ F_T(t)=\Phi\left(\frac{\mu_\beta t-l}{\sqrt{\sigma_\beta^2 t^2+\sigma_B^2 t}}\right)+\exp \left\{\frac{2 \mu_\beta l}{\sigma_\beta^2}+\frac{2 \sigma_B^2 l^2}{\sigma_B^2}\right\} \Phi\left(\frac{2 \sigma_\beta^2 l t+\sigma_B^2\left(\mu_\beta t+l\right)}{\sigma_B^2 \sqrt{\sigma_\beta^2 t^2+\sigma_B^2 t}}\right) \]

\[ \begin{aligned} f_T(t) & =\frac{\alpha^\beta \Gamma\left(\beta+\frac{1}{2}\right)}{\Gamma(\beta)} \frac{l}{\sqrt{2 \pi t^3(1+\eta t)}}\left[\frac{(l-v t)^2}{2 \alpha t(1+\eta t)}+\alpha\right]^{-\beta-\frac{1}{2}} \\ & =\frac{\Gamma\left(\beta+\frac{1}{2}\right)}{\Gamma(\beta)} \frac{l}{\sqrt{2 \pi \alpha t^3(1+\eta t)}}\left[\frac{(l-v t)^2}{2 \alpha t(1+\eta t)}+1\right]^{-\beta-\frac{1}{2}} \end{aligned} \]

由Wiener过程的性质，知道\(\Delta X_{i j}| \mu \sim N\left(\mu \Delta t_{i j}, \sigma^2 \Delta t_{i j}\right)\)，如上面模型类似，似然函数为：

\[ \begin{aligned} l &= \sum_{i=1}^{n} \sum_{j=1}^{m_i} \log \int_{-\infty}^{\infty} f(x;\mu,\sigma) f(\mu) d \mu\\ & = \sum_{i=1}^{n} \sum_{j=1}^{m_i} \log \int_{-\infty}^{\infty}N\left(\Delta x_{i j}; \mu \Delta t_{i j}, \sigma^2 \Delta t_{i j}\right) \times N(\mu;\mu_\beta,\sigma^2_\beta) d\mu \\ \end{aligned} \]

对数似然函数为

\[ \begin{aligned} l &= \sum_{i=1}^{n} \sum_{j=1}^{m_i} \log \int_{-\infty}^{\infty} f(x;\mu,\sigma) f(\mu) f(\sigma) d \mu d\sigma \\ \end{aligned} \]

恒定应力加速Wiener退化过程考虑以下三个假设：

确定正常应力水平\(S_0\)和\(k\)个加速应力水平\(S_1,S_2,\ldots,S_r\)，且满足如下关系 \[{S_00\) 的平稳 Gamma 过程。

产品寿命

因为Gamma 过程的退化是单调递增的，即\(X(t) \sim Ga(\alpha t,\beta)\) ,因此产品寿命的分布可以直接通过退化量的转换获得。当给定产品退化失效阈值 \(\ell\) 时，产品寿命 \(T\) 的CDF和PDF为

\[ \begin{aligned} {F}_{{T}}(t) &={P}({T} \leqslant {t})={P}({X}({t}) \geqslant \ell) \\ &=\int_{\ell}^{\infty} \frac{1}{\Gamma(\alpha {t}) \beta^{\alpha {t}}} {x}^{\alpha {t}-1} \mathrm{e}^{-\frac{{x}}{\beta}} {dx} \\ &=\frac{1}{\Gamma(\alpha {t})} \int_{\frac{\ell}{\beta}}^{\infty} \xi^{\alpha {t}-1} \mathrm{e}^{-\xi} {d} \xi ,\\ {f}_{{T}}(t) &=\frac{{d}}{{dt}} \frac{\Gamma(\alpha {t}, 1 / \beta)}{\Gamma(\alpha {t})} \\ &=\frac{\alpha}{\Gamma(\alpha {t})} \int_0^{\nu / \beta}\left[\ln (\xi)-\frac{\Gamma^{\prime}(\alpha {t})}{\Gamma(\alpha {t})}\right] \xi^{\alpha {t}-1} \mathrm{e}^{-\xi} {d} \xi. \end{aligned} \]

由于\({f}_{{T}}({t})\) 相当复杂，因此一般用B-S分布来逼近：

\[ {F}(t ; \ell)=\Phi\left[\frac{1}{{v}}\left(\sqrt{\frac{{t}}{{u}}}-\sqrt{\frac{{u}}{{t}}}\right)\right], \quad {t}>0. \] 其中，\(\Phi(\cdot)\) 为标准正态分布，\(v=\sqrt{\frac{\beta}{\ell}}, u=\frac{\ell}{\beta \alpha}\). 相应的 PDF 为 \[ {f}(t ; \ell)=\frac{1}{2 \sqrt{2 \pi} {uv}}\left[\left(\frac{{u}}{{t}}\right)^{\frac{1}{2}}+\left(\frac{{u}}{{t}}\right)^{\frac{3}{2}}\right] \exp \left[-\frac{1}{2 {v}^2}\left(\frac{{t}}{{u}}-2+\frac{{u}}{{t}}\right)\right], {t}>0. \]

剩余寿命

根据Gamma 过程的增量独立性，可以得到剩余寿命的分布 \[ \begin{aligned} {F}(t \mid s) &={P}\left({T} \leqslant t \mid {X}(s)={x}_{s}\right) \\ &={P}\left(X(t)-X(s) \geqslant \ell-x_{s}\right) \\ &=\frac{\Gamma\left(\alpha(t-s),\left(\ell-x_{s}\right) / \beta\right)}{\Gamma(\alpha(t-s))} . \end{aligned} \] 这与将失效阈值由 \(\ell\) 变为 \(\ell-{x_s}\) 、时间 \(t\) 变为 \(t-s\) 情形的寿命分布是一样的. 因此，在给定当前状态情况下，可以按照类似的途径更新剩余寿命.

假设共有 \({n}\) 个样品进行性能退化试验. 对于样品 \(i\), 初始时刻 \(t_0\) 性能退化量为 \(X_{i 0}=0\), 在时刻 \(t_1, \cdots, t_{m_i}\) 测量产品的性能退化量, 得到其测量值为 \(X_{i 1}, \cdots, X_{i m_i}\). 记 \(\Delta x_{i j}=X_{i j}-X_{i(j-1)}\) 是样品 \(i\) 在时刻 \(t_{j-1}\) 和 \(t_j\) 之间的性能退化量, 由 Gamma 过程的独立增量性得到\(\Delta X(t) \sim Ga(\alpha \Delta t_{ij},\beta)\)。

由 Gamma 过程独立增量特性, 以及 \[ \Delta X_{i j} \sim G a\left(\alpha \Delta t_{i j}, \beta\right)=\frac{\left(\Delta x_{i j} / \beta\right)^{\alpha \Delta t_{i j}-1}}{\beta \Gamma\left(\alpha \Delta t_{i j}\right)} \mathrm{e}^{-\Delta x_{i j} / \beta} \] 可以获得对数似然函数为 \[ l(\alpha, \beta)=\sum_{i=1}^n\left(\sum_{j=1}^{m_i}\left(\alpha \Delta t_{i j}-1\right) \ln \Delta x_{i j}-\alpha t_{i m_i} \ln \beta-\sum_{j=1}^{m_i} \ln \Gamma\left(\alpha \Delta t_{i j}\right)-\frac{x_{i m_i}}{\beta}\right) \] 由极大似然估计原理, 令 \[ \left\{\begin{array}{l} \frac{\partial l}{\partial \alpha}=\sum_{i=1}^n \sum_{j=1}^{m_i} \Delta t_{i j}\left(\ln x_{i j}-\psi\left(\alpha \Delta t_{i j}\right)-\ln \beta\right)=0 \\ \frac{\partial l}{\partial \beta}=\sum_{i=1}^n\left(\frac{x_{i m_i}}{\beta^2}-\frac{\alpha t_{i m_i}}{\beta}\right)=0 \end{array}\right. \] 式中 \(\psi(u)\) 是 digamma 函数。 \[ \begin{aligned} &\hat{\beta}=\frac{\sum_{i=1}^n X_{i m_i}}{\alpha \sum_{i=1}^n t_{i m_i}}\\ &\sum_{i=1}^n\left[\sum_{j=1}^{m_i} \Delta t_{i j} \ln \left(\Delta x_{i j}\right)-t_{i m_i} \ln \frac{x_{i m_i}}{\alpha t_{i m_i}}-\sum_{j=1}^{m_i} \Delta t_{i j} \psi\left(\alpha \Delta t_{i j}\right)\right]=0 \end{aligned} \]

恒定应力加速Gamma退化过程考虑以下三个假设：

确定正常应力水平\(S_0\)和\(k\)个加速应力水平\(S_1,S_2,\ldots,S_r\)，且满足如下关系 \[{S_0

虽然 Wiener 过程和 Gamma 过程在退化建模中得到了广泛的应用，但这两种模型在实际应用中无法拟合大量退化数据。当这两个过程失败时，另一个选择是逆高斯(inverse Gaussian, IG)过程。

设\(X(t)\)为产品在\(t\)时刻的退化特性。假设退化路径\(X(t)\)可以用参数\(\mu\)和\(\lambda\)的IG过程来建模。此时，\(X(t)\)具有以下属性:

IG 分布的概率密度函数(pdf) \(I G(a, b)\)为 \[ f_{I G}(x, a, b)=\left[\frac{b}{2 \pi x^3}\right]^{1 / 2} \exp \left\{-\frac{b(x-a)^2}{2 a^2 x}\right\} \] 累计分布函数（cdf）为 \[ F_{IG}(x, a, b)=\Phi\left[\sqrt{\frac{b}{x}}\left(\frac{x}{a}-1\right)\right]+\exp \left(\frac{2 b}{a}\right) \Phi\left[-\sqrt{\frac{b}{x}}\left(\frac{x}{a}+1\right)\right], X>0 . \]

产品寿命/首达时分布

令 \(T\) 为\(X(t)\)越过临界值\(l\)的首达时。此时，我们有 \[ T=\inf \{t \mid X(t) \geq l, t>0\} . \] 对于给定的\(l\)，生命周期\(T\)具有以下 cdf

\[ \begin{aligned} F(t) &=P(T \leq t)=P(X(t)>l)=1-F_{I G}\left(l ; \mu t, \lambda t^2\right) \\ &=\Phi\left[ \sqrt{\frac{\lambda}{l}}(t - \frac{l}{\mu}) \right] - \exp\left({\frac{2 \lambda t}{\mu}}\right) \Phi\left[-\sqrt{\frac{\lambda}{l}}(\frac{l}{\mu} +t)\right] \end{aligned}\]

其中 \(\Phi(\cdot)\) 为标准正态分布的cdf。

剩余使用寿命

假设其运行到时刻 \(\tau\) 仍末失效, 且当前性能退化量为 \(x_\tau\left(x_\tau